@keewis,
Instead of changing the code of neighbours_in_kth_ring, would you add a method k_neighbours or kth_neighbours_matrix (or any other better name).
(I plan to test the new method pushing in a vec the result of neighbours_in_kth_ring for k in 0..=k, , and post-sorting it).

Also, I known the difference will be negligible, but to fit with the code style, and since: 1: I don't know the cost of pow(2) and 2: I am not sure the compiler replaces 2 * k by k << 1, would you replace:

(2 * (k  as usize) + 1).pow(2)

by

let k = k as usize;
( (k*(k + 1)) << 2 ) | 1 // = (2k+1)^2

Also, I known the difference will be negligible, but to fit with the code style

sure, I just would have been able to come up with that myself ^^

Instead of changing the code of neighbours_in_kth_ring, would you add a method k_neighbours or kth_neighbours_matrix (or any other better name).

While we're on the subject of naming, would it make sense to call the existing function neighbours_on_kth_ring to make it a bit more apparent that it only returns the pixels on the boundary? Or does the word "ring" already carry that meaning for you?

The new method would be an extension of the neighbours so we could also extend that, with a default of k=1 (but your choice, a new method like neighbours_within_kth_ring, kth_ring_neighbours or, following h3ronpy's terminology, neighbours_disk would be fine with me, as well).

I plan to test the new method pushing in a vec the result of neighbours_in_kth_ring for k in 0..=k, , and post-sorting it

feel free to push to this branch (if I remember correctly, you can do that by adding my fork as a new git remote and fetch it, after which you can check out the branch and push to the remote).

however, if I'm reusing the existing neighbours_in_kth_ring it will take up more memory than the current implementation. Would you mind if I moved the main logic of the method to a helper method? Something like the path_along_cell_side_internal method you already have.

Edit: the recent commits implement this, plus some documentation changes.

Thank you for the implementation @keewis.

Just a few remarks:

If I am correct, the boolean in 1 => self.neighbours(hash, true).values_vec()
should be back to false, and the loop for r in 1..(k + 1) should start at 0.
Also, to avoid an addition, I am pretty sure clippy would say that you should replace ..(k + 1) by ..=k, leading to for r in 0..=k.

You could have kept the first implementation of neighbours_disk (without calling neighbours_in_kth_ring if I remember correctly) and made the current version of neighbours_disk (the one calling neighbours_in_kth_ring) private for testing purposes only: I guess the first version was faster than having to call several times neighbours_in_kth_ring.
But it is ok if you prefer to let as it is now.

I am going to be in holidays for 1 week, I will probably merge the PR when I am back. No hurry?

If I am correct, the boolean in 1 => self.neighbours(hash, true).values_vec()
should be back to false

true, I'll fix that

and the loop for r in 1..(k + 1) should start at 0

No need to call the function for r=0, because that would just return the input cell id. For r=0 the code contains the unconditional result.push(hash) before the loop.

I am pretty sure clippy would say that you should replace ..(k + 1) by ..=k, leading to for r in 0..=k

wow, I didn't know that existed. I agree, that is much more concise.

I guess the first version was faster than having to call several times

If I understand correctly, with neighbours_in_kth_ring_internal there shouldn't be much of a difference now (unless you count the cost of function calls)

I will probably merge the PR when I am back. No hurry?

I've taken the liberty of creating a wrapper package for python (healpix-geo) that depends on the version of cds-healpix-rust of this PR, so we can take our time here. (as pointed out in the OP healpix-geo will be reused for geo-specific things)

(as an additional reason for waiting, we're currently missing a test for neighbours_disk)